Load and analyze results from passive stimulus experiments¶
- Prior to receiving any behavior training, head fixed mice were exposed to three stimulus conditions in randomized order:
- Full screen, square wave gratings, flashed for 200 ms ("visual" condition)
- White noise bursts, 200 ms duration ("auditory" condition)
- Concurrent gratings/white noise bursts ("visual_and_auditory" condition)
- Each experiment was 10-15 minutes in length with ~50 presentations per condition
Imports¶
- The
visual_behavior_researchpackage spun off of the mainvisual_behavior_analysispackage and contains code for processing Inscopix data (https://github.com/AllenInstitute/visual_behavior_research/tree/master/visual_behavior_research) - Helper functions for this analysis are in ..utilities.helper_functions.py
In [167]:
import visual_behavior_research.utilities as vbu
import visual_behavior_research.plotting as vbp
import visual_behavior_research.projects.tbd.database_utilities as du
import visual_behavior_research.projects.tbd.utilities as tbdu
from visual_behavior_research.projects.tbd.session import Session
import imageio
import isx
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import gridspec
import os
import seaborn as sns
# add the top level directory to allow it to be imported
import sys
sys.path.append(os.path.split(os.getcwd())[0])
from utilities import plotting_utilities
from utilities import helper_functions
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
import warnings
warnings.filterwarnings("ignore")
%load_ext autoreload
%autoreload 2
In [2]:
responsiveness_window = 1.5 # number of seconds before/after each stimulus event over which to integrate activity for the purpose of calling a cell 'responsive'
significance_level = 0.05 # level at which to declare a response 'significant'
load summary data¶
In [3]:
summary = tbdu.get_db_table('summary')
mice_table = tbdu.get_db_table('mice')
In [4]:
mice_table
Out[4]:
In [5]:
def get_gender(mouse_id):
mouse_info = vbu.query_labtracks(mouse_id)
return mouse_info['sex']
get sessions marked as 'passive_stim' that are good for analysis¶
In [6]:
all_passive = summary[
(summary['experiment_type'] == 'passive_stim')
&(summary['good_for_analysis']!='False')
# &(summary['anterior_or_posterior']=='anterior')
].copy().sort_values(by=['mouse_id','experiment_type'])
# merge in lens location data from mouse table
all_passive = all_passive.merge(mice_table[['mouse_id','anterior_or_posterior', 'AP_coord', 'ML_coord', 'DV_coord']],on='mouse_id').sort_values(by='recording_date')
# keep only anterior sessions:
all_passive = all_passive[
(all_passive['anterior_or_posterior']=='anterior')
]
# keep only last experiment for each mouse (in some cases we repeated passive stim experiments)
all_passive = all_passive[all_passive[['mouse_id']].duplicated(keep='last')==False].reset_index()
# add a genotype shorthand column for display
all_passive['genotype_shorthand'] = all_passive['genotype'].map(lambda x:helper_functions.get_genotype_shorthand(x))
# sort and reset indices
sessions_to_analyze = all_passive.sort_values(by=['genotype','anterior_or_posterior']).reset_index()
sessions_to_analyze['sex'] = sessions_to_analyze['mouse_id'].map(lambda mid:get_gender(mid.lstrip('M')))
cols_to_display = ['genotype_shorthand','genotype','mouse_id','AP_coord', 'ML_coord', 'DV_coord','sex', 'number_accepted_cells','full_path']
sessions_to_analyze[cols_to_display]
Out[6]:
In [7]:
cols_to_display = ['genotype_shorthand','genotype','mouse_id','AP_coord', 'ML_coord', 'DV_coord','sex']
sessions_to_analyze[cols_to_display].to_csv('/home/dougo/tbd_imaging_mouse_list_2020.09.22.csv', index=False)
In [8]:
sessions_to_analyze['genotype'].value_counts()
Out[8]:
load each mouse's session data into a dictionary called 'session' with mouse_id as key¶
This uses code in visual_behavior_research to build a dedicated session object with many associated methods and attributes
In [9]:
session = {}
for idx,row in sessions_to_analyze.iterrows():
session[row['mouse_id']] = helper_functions.load_session(row['full_path'], load_cached_traces_table=True)
print('on session {} of {}'.format(idx+1,len(sessions_to_analyze)),end='\r')
Apply inclusion criteria to cells in the dataset¶
- keep cells with an SNR above 5 or the 95th percentile of identified peaks in the z-scored data exceeding 3 SDs
- keep cells that have not been labeled as 'likely repeats'
In [10]:
mice = session.keys()
for mouse in mice:
# get metrics dataframe
metrics = session[mouse].metrics
# apply filters, get names of cells to keep
cells_to_keep = metrics[
(metrics['snr']>3)
&(metrics['likely_repeat']==False) # exclude repeats
].reset_index()['cellName']
# save the list of cells to keep as an attribute of the session object
session[mouse].filtered_cell_IDs = cells_to_keep.values
Look at some example cells:¶
a cell with a high SNR¶
In [11]:
mouse = 'M351181'
cell = 'C018'
helper_functions.plot_cell(session[mouse],cell,timeseries_type='filtered')
a cell with a low SNR but a large peaks (the detected events from the Inscopix event detection algorithm don't correspond to the peaks)¶
In [12]:
mouse = 'M432805'
cell = 'C000'
helper_functions.plot_cell(session[mouse],cell)
a cell with low SNR¶
In [13]:
mouse = 'M432805'
cell = 'C002'
helper_functions.plot_cell(session[mouse],cell,timeseries_type='filtered')
However, even an SNR near 3 doesn't guarantee clean calcium transients. It's possible that a higher threshold should be applied to be more conservative.¶
In [14]:
mouse = 'M351181'
cell = 'C150'
helper_functions.plot_cell(session[mouse],cell,timeseries_type='filtered')
I'm also filtering out cells that look to be repeats due to over-segmentation¶
In [15]:
mouse = 'M351181'
session[mouse].repeat_roi_analysis.pairwise_cell_metrics.query('likely_repeat == True')
Out[15]:
In [16]:
mouse = 'M351181'
helper_functions.plot_cell(session[mouse],cell='C079',timeseries_type='raw')
helper_functions.plot_cell(session[mouse],cell='C080',timeseries_type='raw')
helper_functions.plot_cell(session[mouse],cell='C098',timeseries_type='raw')
helper_functions.plot_cell(session[mouse],cell='C097',timeseries_type='raw')
Here is the summary of sessions and cell counts before/after filtering (sessions with incomplete data aren't included)¶
In [17]:
sessions_to_analyze
Out[17]:
In [18]:
sessions_to_analyze['cell_count_post_filtering'] = sessions_to_analyze['mouse_id'].map(lambda x:len(session[x].filtered_cell_IDs))
sessions_to_analyze.set_index('mouse_id',inplace=True)
cols_to_display = ['genotype_shorthand','anterior_or_posterior','number_accepted_cells','cell_count_post_filtering']
sessions_to_analyze[cols_to_display]
Out[18]:
for each session, calculate responsiveness data for each accepted cell¶
extract the integral of the z-scored response in the
In [19]:
mice
Out[19]:
In [20]:
for mouse_id in mice:
session[mouse_id].responsiveness_data = helper_functions.get_responsiveness_data(session[mouse_id], window_size=responsiveness_window, behavior_condition='passive')
In [21]:
mouse_id = 'M335379'
session[mouse_id].responsiveness_data.sample(15,random_state=0).set_index('cell_id')
Out[21]:
for each session, build responsiveness summary for each accepted cell¶
For each cell:
- build a distribution of post-stimulus response magnitudes for a given stimulus condition
- build a distribution of pre-stimulus response magnitutes for every stimulus condition
- test the hypothesis that these distributions are different (i.e., a cells calcium activity in the 1 second following the stimulus differs from its activity in the spontaneous state) using an independent t-test
- save the p-value for each comparison
In [22]:
for mouse_id in mice:
session[mouse_id].responsiveness_summary = helper_functions.build_responsiveness_summary(session[mouse_id])
In [23]:
mouse_id = 'M375765'
session[mouse_id].responsiveness_summary.sort_values(by='p-value_on_integrals').head(10)
Out[23]:
In [24]:
mouse_id = 'M375765'
helper_functions.plot_cells_with_significant_responses(session[mouse_id]);
Here's an example of how the significance calculation works:¶
In [25]:
fig=plt.figure(figsize=(4,3))
ax = []
ax.append(vbp.placeAxesOnGrid(xspan = (0,1),yspan=(0,0.5),fig=fig))
ax.append(vbp.placeAxesOnGrid(xspan = (0.25,0.75),yspan=(0.6,1),fig=fig))
window_size = responsiveness_window
passive_stim_trials = session['M334410'].behavior_core_data['trials']
stim_condition = 'auditory'
dat = helper_functions.get_responses(
session['M334410'],
'C071',
passive_stim_trials[passive_stim_trials['stim_type'] == stim_condition]['nearest_F_frame']
)
pre_stim_indices = np.where(np.logical_and(dat['t']>=-window_size,dat['t']<0))
post_stim_indices = np.where(np.logical_and(dat['t']>0,dat['t']<=window_size))
pre_stims = []
post_stims = []
for idx in range(len(dat['all_traces'])):
pre_stims.append(np.mean(dat['all_traces'][idx][pre_stim_indices]))
post_stims.append(np.mean(dat['all_traces'][idx][post_stim_indices]))
df = pd.DataFrame({
'value':pre_stims+post_stims,
'condition':['pre_stim']*len(pre_stims) + ['post_stim']*len(post_stims)
})
vbp.plot_event_triggered_timeseries(dat,ax=ax[0])
ax[0].axvspan(-window_size,0,color='orange',alpha=0.25)
ax[0].axvspan(0,window_size,color='blue',alpha=0.25,zorder=-1)
ax[0].set_xlim(-2.5,2.5)
c=sns.stripplot(
x='condition',
y='value',
data=df,
ax=ax[1],
palette=sns.color_palette(['orange', 'blue']),
alpha=0.5
)
b=sns.boxplot(
x='condition',
y='value',
data=df,
ax=ax[1],
palette=sns.color_palette(['gray','gray']),
whis=np.inf
);
ax[0].set_ylim(-2,6)
ax[1].set_ylim(-2,6)
sns.despine()
fig.tight_layout()
How dependent is the significance value on the size of the window?¶
Look at four cells for mouse 375765
In [26]:
mouse_id = 'M375765'
for cell_id in ['C009','C010','C012','C006']:
fig,ax=helper_functions.plot_significance_vs_window_size(mouse_id,cell_id,session,metric='mean')
what is the overall fraction of responsive cells?¶
In [27]:
for mouse_id in mice:
session[mouse_id].responsive_cells = {}
responsive_cells = helper_functions.get_cells_with_significant_responses(session[mouse_id], significance_level=significance_level)
sessions_to_analyze.at[mouse_id,'number_of_responsive_cells'] = len(responsive_cells)
passive_stim_trials = session[mouse_id].behavior_core_data['trials']
for cell_id in responsive_cells:
session[mouse_id].responsive_cells[cell_id] = {}
for col,condition in enumerate(['visual','auditory','auditory_and_visual']):
dat = helper_functions.get_responses(
session[mouse_id],
cell_id,
passive_stim_trials[passive_stim_trials['stim_type'] == condition]['nearest_F_frame']
)
d = {'trial_{:0>2d}'.format(i):arr for i,arr in enumerate(dat['all_traces'])}
d['t'] = dat['t']
session[mouse_id].responsive_cells[cell_id][condition] = pd.DataFrame(d).set_index('t')
def get_fraction_responsive(row):
return row['number_of_responsive_cells']/row['cell_count_post_filtering']
sessions_to_analyze['fraction_responsive_cells'] = sessions_to_analyze.apply(get_fraction_responsive, axis=1)
In [28]:
cols_to_display = ['genotype_shorthand','anterior_or_posterior','number_accepted_cells','cell_count_post_filtering','number_of_responsive_cells','fraction_responsive_cells']
sessions_to_analyze[cols_to_display]
Out[28]:
make some summary plots¶
lump cells together, make an overall summary plot¶
In [29]:
mice
Out[29]:
In [30]:
responsiveness_summary_all = [pd.DataFrame()]
responsiveness_data_all = [pd.DataFrame()]
for mouse_id in mice:
mouse_summary = session[mouse_id].responsiveness_summary
mouse_summary['mouse_id'] = mouse_id
mouse_summary['location'] = sessions_to_analyze.query('mouse_id == @mouse_id')['anterior_or_posterior'].item()
mouse_data = session[mouse_id].responsiveness_data
mouse_data['mouse_id'] = mouse_id
mouse_data['location'] = sessions_to_analyze.query('mouse_id == @mouse_id')['anterior_or_posterior'].item()
responsiveness_summary_all.append(mouse_summary)
responsiveness_data_all.append(mouse_data)
responsiveness_summary_all = pd.concat(responsiveness_summary_all)
responsiveness_data_all = pd.concat(responsiveness_data_all)
def make_unique_cell_id(row):
return row['mouse_id'] + '_' + row['cell_id']
responsiveness_summary_all['unique_cell_id'] = responsiveness_summary_all.apply(make_unique_cell_id,axis=1)
In [31]:
responsiveness_summary_all['integral_change_magnitude'] = responsiveness_summary_all['post_stim_integral_mean'] - responsiveness_summary_all['pre_stim_integral_mean']
responsiveness_summary_all['mean_change_magnitude'] = responsiveness_summary_all['post_stim_mean_all'] - responsiveness_summary_all['pre_stim_mean_all']
In [32]:
responsiveness_summary_all.head()
Out[32]:
In [33]:
responsiveness_summary_all['mean_change_magnitude'].describe()
Out[33]:
In [34]:
def plot_mod_index_hist(index_vals, significant, ax, bin_edges=np.arange(-1,1,0.05)):
'''
index_vals (array of float): values of the index you want to plot
significant (array of bool): whether each index val is significant or not
'''
ax.hist([index_vals[significant],index_vals[~significant]],
color=['k','darkgrey'], edgecolor='None', stacked=True,
bins=bin_edges, rwidth=0.8, normed=True)
def plot_response_histograms(ax=None, significance_threshold=0.05):
'''
ax should be a 1x3 array of axis handles if passed
'''
if ax is None:
fig,ax=plt.subplots(1,3,figsize=(12,2),sharex=True, sharey=True)
for ii,condition in enumerate(['visual','auditory','auditory_and_visual']):
magnitude = responsiveness_summary_all.query('condition == @condition')['mean_change_magnitude']
significant = responsiveness_summary_all.query('condition == @condition')['p-value_on_means'] < significance_threshold
plot_mod_index_hist(
magnitude,
significant,
ax[ii]
)
axins = inset_axes(ax[ii], width="40%", height="40%", loc=1)
axins.pie(
[significant.sum(),len(significant)-significant.sum()],
colors=['black','gray'],
labels=['','']
)
ax[ii].set_title('condition = {}\nfraction modulated cells: {:0.3f}'.format(condition,significant.mean()))
ax[ii].set_xlabel('change in mean (z-score)', fontsize=16)
ax[ii].set_xticks([-1, -0.5, 0, 0.5, 1])
ax[ii].set_yticks([])
ax[0].set_ylabel('fraction')
sns.despine()
In [35]:
plot_response_histograms()
In [163]:
def plot_all(mouse_id,cell_id,ax=None, label_x=True, show_individual_traces=False):
sess = session[mouse_id]
if ax is None:
fig,ax = plt.subplots(1,3,figsize=(8,4),sharex=True,sharey=True)
titles = [
'visual stimuli\np={:0.5f}',
'auditory stimuli\np={:0.5f}',
'visual+auditory stimuli\np={:0.5f}',
]
rs = responsiveness_summary_all.query("mouse_id==@mouse_id and cell_id==@cell_id")[['condition','p-value_on_means']]
passive_stim_trials = sess.behavior_core_data['trials']
for col,cond in enumerate(['visual','auditory','auditory_and_visual']):
dat = helper_functions.get_responses(
session[mouse_id],
cell_id,
passive_stim_trials[passive_stim_trials['stim_type'] == cond]['nearest_F_frame']
)
response_df = []
for ii in range(np.shape(dat['all_traces'])[0]):
df = pd.DataFrame({'t':dat['t'],'z_scored_activity':dat['all_traces'][ii]})
df['condition'] = cond
df['repeat_number'] = ii
response_df.append(df)
response_df = pd.concat(response_df)
if show_individual_traces:
ax[col] = vbp.plot_event_triggered_timeseries(
dat,
ax=ax[col],
foreground_color='black',
)
else:
sns.lineplot(
data = response_df,
x='t',
y='z_scored_activity',
ax = ax[col],
color='black',
)
ax[col].axvline(0,color='k',linewidth=3,zorder=-np.inf,alpha=0.5)
ax[col].set_xlim(-1.5,3)
ax[col].set_title(titles[col].format(rs[rs['condition']==cond]['p-value_on_means'].item()))
if col == 0:
ax[col].set_ylabel('z-scored activity',fontsize=14)
ax[col].set_yticklabels(ax[col].get_yticks(),fontsize=12)
if label_x == False:
ax[col].set_xticks([])
ax[col].set_xticklabels([])
ax[col].set_xlabel('')
else:
ax[col].set_xlabel('time from event (s)')
sns.despine()
In [165]:
mouse = 'M375765'
cell_id = 'C009'
plot_all(mouse ,cell_id, label_x=True)
In [38]:
mouse = 'M351181'
cell_id = 'C104'
plot_all(mouse,cell_id)
In [39]:
responsiveness_summary_all[
(responsiveness_summary_all['condition']=='auditory')
&(responsiveness_summary_all['p-value_on_means']<significance_level)
]
Out[39]:
In [40]:
def show_schematic(ax=None):
if ax is None:
fig,ax=plt.subplots()
image = imageio.imread('passive_stim_schematic.png')
# trim off black border on PNG
ax.imshow(image[10:-10,10:-10,:],interpolation='nearest')
ax.axis('off')
In [41]:
fig,ax=plt.subplots()
show_schematic(ax)
Final Figure¶
In [42]:
%matplotlib inline
In [166]:
save_fig = True
figsize=(2*8.5,2*6)
figure_2 = plt.figure(figsize=figsize)
axes = {
'A':plotting_utilities.placeAxesOnGrid(figure_2, xspan=[0.05,0.4], yspan=[0,0.6]),
'B':plotting_utilities.placeAxesOnGrid(figure_2, xspan=[0.5,1], yspan=[0.1,0.28], dim=[1,3], sharex=True, sharey=True, wspace=0.02),
'C':plotting_utilities.placeAxesOnGrid(figure_2, xspan=[0.5,1], yspan=[0.34,0.52], dim=[1,3], sharex=True, sharey=True, wspace=0.02),
'D':plotting_utilities.placeAxesOnGrid(figure_2, xspan=[0.05,1], yspan=[0.67,1], dim=[1,3], sharex=True, sharey=True, wspace=0.02),
}
show_schematic(axes['A'])
plot_all('M375765', cell_id='C009', ax=axes['B'], label_x=False)
plot_all('M351181', cell_id='C104', ax=axes['C'], label_x=True)
plot_response_histograms(axes['D'])
# add labels
labels = [
{'label_text':'A', 'xpos':0, 'ypos':0.05},
{'label_text':'B', 'xpos':0.4, 'ypos':0.05},
{'label_text':'C', 'xpos':0.4, 'ypos':0.30},
{'label_text':'D', 'xpos':0, 'ypos':0.60},
]
for label in labels:
plotting_utilities.add_label(
figure_2,
label_text = label['label_text'],
xpos = label['xpos'],
ypos = label['ypos'],
fontsize=20,
weight="bold"
)
if save_fig:
plotting_utilities.save_figure(figure_2, fname='figure_2', formats=['.png','.pdf'], size=figsize, dpi=500)
heatmaps of individual cells¶
In [157]:
def make_heatmaps(mouse_id, cell_id):
fig,ax = plt.subplots(1,3,figsize=(18,4))
for col,stim_condition in enumerate(['visual','auditory','auditory_and_visual']):
p_val = responsiveness_summary_all.query('mouse_id == @mouse_id and cell_id == @cell_id and condition == @stim_condition')['p-value_on_means'].iloc[0]
title = 'condition = {}, p-value = {:0.4f}'.format(stim_condition, p_val)
helper_functions.single_cell_heatmap(
session = session[mouse_id],
cell_id = cell_id,
events = session[mouse_id].behavior_core_data['trials'].query('stim_type == @stim_condition')['nearest_F_frame'],
ax=ax[col],
cbar = True,
title=title
)
fig.tight_layout()
plt.subplots_adjust(top=0.85)
fig.suptitle('mouse_id = {}, cell_id = {}'.format(mouse_id, cell_id),fontweight='bold')
In [159]:
for idx,row in responsiveness_summary_all[responsiveness_summary_all['p-value_on_means'] < significance_level].drop_duplicates('unique_cell_id').sort_values(by='p-value_on_means').iterrows():
make_heatmaps(mouse_id = row['mouse_id'], cell_id = row['cell_id'])
In [ ]: